Single Molecule Chemical Reaction: Kramers Approach Revisited 



in 
o 
o 

(N 
5h . 

< 
OO 



CO 



I 

C 

O 

o 



> 
in 
in 
^t- 
^l- 
o 
in 
o 



-a 
c 

o 
o 



X 
S3 



G. Margolin 1 and E. Barkai 1 ' 2 

'Department of Chemistry and Biochemistry, Notre Dame University, Notre Dame, IN 46556 
2 Department of Physics, Bar Ran University, Ramat Gan, Israel 52900 
(Dated: 2nd February 2008) 

Single molecule chemical reactions yield new insight into fluctuation phenomena which are obscured in mea- 
surement of ensemble of molecules. Kramers escape problem is investigated here in a framework suitable for 
single molecule reactions. In particular we obtain distributions of escape times in simple limiting cases, rather 
than their mean, and investigate their sensitivity on initial conditions. Rich physical behaviors are observed: sub- 
Poissonian statistics when the dynamics is only slightly deviating from Newtonian, super-Poissonian behavior 
when diffusion is dominating, and Poissonian behavior when Kramers original conditions hold. By varying 
initial conditions escape time distributions can follow a (usual) exponential or a r~ 3//2 decay, due to regular 
diffusion. We briefly address experimental results which yield the r _3//2 behavior (with cutoffs) and propose 
that this behavior is universal. 
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Chemical reaction of a single molecule evolving between 
two states A ^ B or of two species A + B ^ AB is now 
followed in many laboratories using single molecule spec- 
troscopy techniques mil Hi SHU- Such experiments 
yield detailed statistical information on chemical conforma- 
tional changes, and simple chemical processes in condensed 
phase environments, for example the distribution of occupa- 
tion times of states A and B in the process A ^ B. Such infor- 
mation is impossible to obtain when measurements are made 
on many molecules, since the ensemble averaging wipes out 
the detailed information found on the single molecule level. 
For ensemble of molecules, usually simple reactions are as- 
sumed to follow a rate process. 

For ensemble of molecules undergoing a chemical reaction, 
classical concepts like reaction coordinate, and rate equations, 
work well in many cases. In particular Kramers model 
for activation over a barrier, is a fundamental tool for model- 
ing chemical reactions in condensed phase environments. In 
this manuscript we analyze Kramers problem, in a framework 
of single molecules. We first discuss Kramers original ap- 
proach, and its limitations in the single molecule domain. 

Kramers describes a chemical reaction in terms of a re- 
action coordinate x(t). The complicated interaction of the 
chemical species with their environment is replaced with a 
stochastic one dimensional approach. The coordinate x(t) 
evolves in a deterministic force field V(x), and is also cou- 
pled to a thermal heat bath. The reaction coordinate is sup- 
posed to escape a metastable state. The inverse of the average 
time of escape (r), from the bottom of the well, serves as an 
estimate on the ensemble averaged reaction rate. Two impor- 
tant regimes of Kramers problem are the underdamped and 
overdamped limits. Many refinements, non-trivial results and 
generalizations of Kramers problem are known, e.g. Kramers 
turn over behavior, quantum effects, and non-Markovian gen- 
eralizations (e.g., IKlll 1 1 1 1 2fl : see |9[ for a review). Experi- 
mental validation of the theory is also obtained fTjll . 

At least three aspects of Kramers problem must be revised 
in the context of single molecule chemical reaction. The most 



obvious one is that now we must consider the distribution of 
occupation times in a chemical state, and not limit the theory 
to averaged escape rates. Previous work considered temporal 
dependence of the rate 1 1411 1511 until it reaches an equilibrium 
value in a single escape (transition) event, or the transient be- 
havior after the particle injection close to the bottom of the 
well 1 16]. The idea of fluctuating rates in multiple transi- 
tions, in the context of single molecule experiments has been 
scrutinized and used in, e.g., Ill Il7l ll 81 Il9ll . A second is- 
sue is the sensitivity of single molecule chemical reactions 
to initial conditions. Consider the ongoing chemical reaction 
A B. For example using fluorescence resonance energy 
transfer (FRET) methods, one may follow the closin g (s tate A) 
and opening (state B) of a large molecule flfsl 161 12(112 ill . The 
experimental data then yields the string of occupation times 



.}. Following Kramers assume that such events 



are described by a reaction coordinate, which goes over a sin- 
gle barrier, to cross from state A — > B and vice versa. Also 
we assume that a measurement may distinguish between state 
A and state B, as is shown in many experiments. In terms 
of the reaction coordinate this means that when x(t) < Xb 
the system is in the state A, otherwise it is in the state B and 
X), is the boundary (usually and conveniently assumed to be 
at the top of the potential barrier). Then immediately after the 
transition event from say A — > B, the reaction coordinate is in 
close vicinity to Xb- Thus a short time after the transition from 
state A — > B there is an increased likelihood for a back transi- 
tion B — > A. This possibility is expected to yield bunching of 
chemical activity on the time axis, i.e. to intermittency where 
strong activity is observed for some period of time followed 
by periods of lesser activity. Of course to observe such effects 
the resolution time in the experiments must be short compared 
to relaxation times of the dynamics. A third important point 
is that in single molecule experiments, an additional length 
scale A is introduced into the problem. For example the radius 
of the laser illumination field, or the Forster length scale in 
FRET experiments |4, 21]. Roughly, in the reaction A ^ B, 
a fluorescence signal is recorded if an acceptor and donor are 
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Figure 1: Escape from a metastable well. In the underdamped case, 
the relevant coordinate is the action I, which is related to the en- 
ergy E of the particle by / = 2nE/u. In the overdamped case, the 
relevant coordinate is position x in real space. 



within a range of Forster radius from each other (state A) while 
otherwise the signal is zero (state B). In terms of the reaction 
coordinate this implies that when x(t) < A the system is in 
the state A otherwise it is in the state B. The point to notice 
is that in principle A can be anywhere along the reaction coor- 
dinate. In particular the classical approach of an escape from 
a bottom of the well over a maximum in the potential field is 
not expected to be general. 

Hence we investigate Kramers escape problem, obtaining 
distribution of escape times. With initial conditions both in 
vicinity of the escape point and far from it. We classify the 
deviations from exponential behaviors, and show that in many 
cases rate concept is not valid. In particular we classify a 
turnover behavior from power law behavior to exponential, 
and show that the power law behavior is general. 

Model 1 We consider a classical particle undergoing un- 
derdamped diffusion in a harmonic potential well. As well 
known the relevant coordinate for underdamped motion is the 
action /. The probability density function (PDF) P of finding 
the particle with action / obeys 



dP 
~dt 





Ti 



ip 



2irk B T dP 
~u ~dl 



(1) 



where lu is the harmonic frequency, E is the particle energy, 
ksT is the thermal energy, and 7 is the damping coefficient; 
7 <C co (see Fig. [J] J = 2ttE/lu for constant u>). In 
what follows we use dimensionless time t = -ft and action 
T = Ioj /ZirksT . Let <f)(r) be the PDF of escape times from 
Zo to T esc > To- Mathematically we use absorbing boundary 
conditions at T esc , so that 4>(t) is the first passage time dis- 
tribution (FPTD) from To to T esc . We obtained the Laplace 
t — > u transform of <P(t), using known solution to Kummer's 
equation We find 



iF\(u; 1;2q) 



(2) 



where iFi(a;b;c) is the regular confluent hypergeometric 
function. Since i.Fi(0; b; c) = 1 the function 4>(t) is normal- 
ized to 1 . We note that a second presentation of the solution 
in terms of an infinite sum of exponentially decaying modes is 
possible, in time domain. However, as we shall show now for 
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Figure 2: Q + 1 as a function of r = To/T esc . Full line is the 
asymptotic line for small Te.sc, given by (1 + r)/(2(l — r)). Dashed 
line is for T ESC — 1 and dot-dashed line is for I e3C — 10. Notice 
the logarithmic scale. Two generic behaviors are observed: (i) if 
To < 1» a smooth transition from sub-Poissonian to super Poissonian 
behavior is observed (ii) if To > 1, a Poissonian behavior is found 
until To — Tesc and then a sharp transition from Poissonian to super- 
Poissonian behavior is observed. 



an important parameter regime, one cannot replace such stan- 
dard eigenfunction solutions with a summation over a finite 
number of modes consisting of a few exponential functions, 
rather the solution exhibits a power law behavior. 

To quantify deviations from exponential statistics, we use 
the parameter 



1. 



where (r) is the average escape time and a 2 = (r 2 ) — (r) 2 
is its variance. If we have Poissonian behavior consistent with 
a rate equation approach, then Q = 0. When the dynamics 
is Newtonian (i.e. diffusion is weak) the PDF of escape time 
is narrow, leading to Q < 0, a sub-Poissonian behavior. On 
the contrary, if the PDF of escape times is wide spread, there 
is a super-Poissonian behavior and Q > 0. Using the small u 
expansion of the exact solution, we find 



(r) = A(I esc ) - A(I ) 



(3) 



and 



with 



a 2 = A 2 (I esc ) - A 2 (I ) - 2[B(I esc ) - B(l )} (4) 
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(5) 



In Fig. |2 we plot Q + 1 versus Tq for three cases where 
T esc < 1 (escape over a shallow barrier), I esc = 1 and 
l esc > 1 (escape over a large barrier). As To — > T esc , 
namely the case when the chemical reaction starts close to 
the escaping zone, we observe Q — > 00, i.e. strong super- 
Poissonian behavior. For T esc < 1 we do not expect, and 
indeed do not find an exponential behavior, since the barrier 
is not high. Here, depending on initial position, either sub- 
or super-Poissonian behavior is generally observed. As To 
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Figure 3: Distribution of escape times 4>( T ) m the underdamped case. 
At early and intermediate times 4>{r) is governed by Eq. ^8), which 
exhibits a power law decay r _3//2 at intermediate times (thin line), if 
the initial condition Xq is close to the escape condition T esc . At later 
times, there is an exponential decay with rate given by Kramers rate. 



becomes closer to T esc , diffusion in Eq. Q becomes more 
dominant than the deterministic drift, and hence Q grows. 
Conversely, when To <§; T esc and T esc < 1 drift becomes 
more important, leading to negative Q. However, diffusion 
can never be neglected as it is the only mechanism leading to 
an eventual escape from the well. Note that always Q > — 1/2 
and the absolute minimum of Q = —1 is unachievable in 
this model. For T esc > 1 and To sufficiently below T esc , we 
have exponential behavior and Q ~ 0. The latter case corre- 
sponds to Kramers' original treatment of escape from a deep 
metastable state. 

We now consider 4>(t). In Laplace space and in the limit 
u ^> max{l, 1/Tq} we use lEill to reduce Eq. to 



Model 2 We turn now to the other limit of particle diffusion 
- namely to the overdamped limit, and demonstrate the devi- 
ation from Poissonian behavior. The relevant coordinate here 
is the spatial coordinate x. Restricting ourselves to one spatial 
dimension, the following equation is then obtained: 



dP(x,t) 
dt 



P(x,t), (9) 



where m is the particle mass and F(x) = —U'(x) is the deter- 
ministic part of the force, due to potential field. We consider 
two cases, when U{x) — ±mu 2 x 2 /2 near the escape point, 
i.e., parabolic and inverted parabolic potential, with to <C 7. 
We define dimensionless time r = t/8 and position y = x/l, 
where — ^uj~ 2 and I = \JksT j (mcj 2 ). 

Parabolic potential. We are interested in the FPTD from 
initial position yo to y esc . This function is known in Laplace 
space and is given by (e.g., |23]) 



0(«) 



D_ u {-yos)eyH A H_ u (-y s /s/2) 



D^ u (-y esc s) e y 2 ^/' 1 H- a (- 



(10) 



where D\(z) is the parabolic cylinder, or Weber function, 
H- U (z) is the generalization of Hermite polynomials used in 
MATHEMATICA and s = sign(z/ esc - y ). The PDF 0(r) is 
normalized to 1. Using various formulas from I22I1 we can 
simplify Eq. (TiPl for large u ^> ma,x{y 2 sc , y 2 }: 



4>{u) ~ exp 



Vesc 



-Vu\y esc -y \ 



^esc 

2rT 



1/4 



= (X -X eS c)/2 -2V^(v / X^T-v / 2o) 



(6) 



The important thing to notice is the non-analytical behavior of 
this asymptotic solution, namely the term of the form e~^". 
This term is an indication for power law behavior (with some 
cutoffs) since if <p{u) = for all u then 4>(t) ~ r~ 3 / 2 for 
large r (e.g., 12411 ). To find r for which (f>(r) ~ r~ 3 / 2 we have 
to demand the validity of Eq. l|6) for small ^/u{\/I esc — \/Tq), 
while u is large. This yields 



< r < min{l,2o} 



(7) 



so that the reaction should start in the vicinity of the escaping 
point. This behavior is demonstrated in Fig. [3] 

When r^Owe find 0(r) — > 0, a non-exponential behav- 
ior. This behavior is due to the fact that it takes the particle a 
certain amount of time to reach the boundary. More quantita- 
tively, for short times Eq. (|6j is always a good approximation 
for (j>(u), leading to 



(r) 



7 

-t--esc 

~2b~ 



1/4 



e (Xo-X csc )/2 (v ^- 



exp[-(VT ese - Vlo) 2 /r] 



(8) 



exhibiting a non-analytical behavior of the e _v/ " type similar 
to the previous example. In order to have the <j)(r) oc r~ 3 / 2 
scaling we thus have to demand the validity of this approxima- 
tion for small y/u\y esc — yo\ < 1. Together with the condition 

u ^> max{y 2 sc , y 2 } it leads to 



{yesc - yo) < T < 



(11) 



As we are interested in the escape over the barrier, we have 
either y Q < y esc if y esc > 0, or y > y esc if yesc < 0. In a 
deep well (low T), if \yo\ -C \y e sc\ men </>(x) will be nearly 
exponential, similar to Model 1. If the well is not deep, or 
yo and y esc are close, Q shows a super Poissonian behavior. 
Noticeably negative Q is only possible if the initial position 
of the particle is on the other side of the well from the escape 
point, and high above it in energy (when yo/yesc becomes 
very negative). 

Inverted parabolic potential. The FPTD is now given by 



4>{u) 



-u(-yos) 



-D-i-u(-y e 

The normalization is 



exp 



Vesc 



(12) 



*«)) = i+erf( ^; < 1. 

1 +erf(y esc s/\/2) 
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The normalization is less than 1 in this case because some 
particles will escape in the direction of — soo, where the po- 
tential drops unbounded. Of course, if the particles are not 
allowed to escape to infinity, so that the only exit is through 
y esc then the normalization should be 1 . However, the above 
formulae serve as a good approximation for a deep well, if 
we consider the behavior around the parabolically shaped es- 
cape barrier at sufficiently short times. Similar to the previous 
section, for u 3> max{l, y\ scl y 2 } 



4>{u) ~ exp 



2 2 
Vesc ~ HO 



V«|t/eae-i/o| 



and the condition for observing r 3 / 2 scaling is 



■Vo) <r< 



max{l,y 2 sc ,y 2 } 



In this Letter, we demonstrated non-Poisson statistics of es- 
cape from a potential well, which contradicts reaction rate ap- 
proach used for ensemble (bulk) dynamics. Both sub- and 
super-Poisson statistics were observed, indicated by the sign 
of the Q parameter. For small travel distances, and for a shal- 
low well in the overdamped case, scaling 4> cx t~ 3 / 2 is dom- 
inant and gives rise to a strong super-Poissonian statistics (cf. 
Fig. |3- The r~ 3 / 2 behavior is a result of diffusion processes 
12511 . The occurrence of such a behavior is easily understood 
if one realizes that in many cases a Fokker-Planck equation 
with variable coefficients (e.g., Eq. {I}) can locally be ap- 
proximated by an equation with constant coefficients, i.e., by 
advection-diffusion equation. As travel distance increases, 
characteristic time grows and drift term can compete with dif- 
fusion, allowing for Poissonian and sub-Poissonian behavior. 
In the underdamped case sub-Poissonian behavior is found in 
many cases (see Fig. |2ji, in the overdamped case we may ob- 
serve sub-Poissonian behavior only under special conditions. 
As is well known, near-exponential 4>{t) is obtained for a deep 
well if initial energy is sufficiently below the escape energy. 

Finally let us compare our model results to experiments on 
single molecules and to other models. Blinking nanocrystals 
exhibit a behavior of occupation times of on and off times 
close to r~ 3 / 2 with cutoffs 1261 1271 12811 . The blinking is be- 
lieved to describe charging of a single nanocrystal (charged 
NC can be off). Shimizu et al. |28] briefly suggested a diffu- 
sion in energy space to describe such behavior, however did 
not consider the effect of dissipation and temperature which 
are always present in case of diffusion in energy space. Re- 
cent experiments on diffusing beads that come in and out of 
focus of a laser field, also exhibit the r -3 / 2 behavior 12^1 . 
Occupation times in single molecule Raman experiments are 
described also by r~ 3 / 2 behavior 1 30] . Dynamics of single 
ion channels sometimes exhibit the r~ 3 / 2 13 ll 13211 . Goychuk 
and Hanggi 1 32] suggested a model based on a reaction where 
space is divided into two: a zone with free diffusion (which 
yields the 3/2 law) and a zone where the reaction coordinate 
is climbing over a potential field. While all these systems and 
models are very different, they all exhibit a universal tendency 



for a t -3 / 2 and in several cases an exponential cutoff is ob- 
served. Thus the turn over behavior we found, from a power 
law behavior, to an exponential behavior, will be a useful con- 
cept for the single molecule domain. Further we do not ex- 
pect this behavior to be limited to Kramers problem, since 
diffusion is expected to control short time dynamics of many 
reactions. 
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